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Abstract 

The correct computation of orbits of discrete dynamical systems on the inter- 
val is considered. Therefore, an arbitrary-precision floating-point approach 
based on automatic error analysis is chosen and a general algorithm is pre- 
sented. The correctness of the algorithm is shown and the computational 
complexity is analyzed. There are two main results. First, the computational 
complexity measure considered here is related to the Lyapunov exponent of 
the dynamical system under consideration. Second, the presented algorithm 
is optimal with regard to that complexity measure. 
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1. Introduction 

Consider a discrete dynamical system {D, /) on some compact interval 
D CW, called the phase space, given by a function f : D D, a recursion 
relation Xn+i = f{xn) and an initial value xq G D. The sequence of 
iterates is called the orbit of the dynamical system in phase space corre- 
sponding to the initial value xq. If such a dynamical system is implemented, 
that is a computer program is written for calculating a finite initial segment 
of the orbit for given xq, care has to be taken in choosing the appropri- 
ate data structure for representing real numbers. Traditionally, IEEE 754 



double floating-point numbers [IJ, |9| are used. However, if the dynamical 
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system shows chaotic behavior, a problem arises. The finite and constant 
length of the significand of a double variable causes rounding errors which 
are magnified after each iteration step. Already after a few iterations, the 
error is so big that the computed values are actually useless. For example in 
(itI [21I this phenomenon is examined for the dynamical system (D, f) with 
D = [0, 1], /(x) = 3.75 ■ X ■ (1 — x) and the initial value Xq = 0.5. To put 
things right, a rigorous method for computations with real numbers has to 
be used. There already exist some rigorous numerical methods in the field 



of dynamical systems and chaos [ll|, [lOl, |30|, l23|, l25|, |20 



In the next section, a rigorous method based on arbitrary-precision fioat- 
ing point arithmetic is presented and used to investigate the iteration of a 
generalization of the above mentioned function. Correctness of the results 
are obtained by using a method called running error analysis. The method 
and the numerics are compared to interval arithmetic. In Section |3l the algo- 
rithm is generalized to arbitrary functions /. The aim of the present paper 
is to give bounds on some kind of space complexity of the algorithm. To be 
more precise, the behavior of the the length of the significand in arbitrary- 
precision arithmetic is analyzed in the task of iterates of discrete dynamical 
systems. The minimal length of the significand needed for fioating-point 
numbers such that any computed point of an initial segment of the orbit has 
a specified and guaranteed accuracy is examined. This minimal length will 
be related to the length of the initial segment of the orbit. To cope with 
this task, a precise mathematical framework for fioating-point computations 
is applied. This framework should be suited to computability concepts over 
the reals. Finally, a complexity measure for describing the computational 
effort on arbitrary-precision fioating-point numbers is introduced. Roughly 
speaking, it is the ratio of the length of the significand to the number of it- 
erations in the limit of number of iterations to infinity. The first main result 
shows that this complexity measure is related to the Lyapunov exponent. 
The second main result proves that the presented algorithm to compute the 
orbit up to any given accuracy is optimal with respect to that complexity 
measure. As a consequence, these results give some advice for economically 
designing reliable algorithms simulating one-dimensional discrete dynamical 
systems. 
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2. Dynamic behavior of the logistic equation and rounding error 



In this section, the discrete dynamical system {D, f^) with D = [0, 1] and 
: D ^ D, f^{x) := — x) for some control parameter /x G (0,4] is 

investigated. In the literature, the recursion relation Xn+i = /^(a^n) is called 
the logistic equation When implementing the logistic equation on a real 
computer and demanding to obtain true values for the orbit {xn)n, some 
rigorous method is needed. Since for some values of /i the dynamics is highly 



chaotic, inaccuracies are magnified exponentially in time [6|, |13|. Therefore, 
it is clear that using floating-point numbers with a predefined, fixed precision 
makes sense only if the maximum iteration time N also is a predefined, fixed 
number. If the algorithm should work for any N, a more elaborate approach 
is needed. First one can work with arbitrarily high precision fioating-point 
numbers, the precision dynamically set and the error control implemented 
in the algorithm separately. A software package for doing this task is for 
example MPFR [sj. Second, there are methods with automatic error control. 



for example interval arithmetic |19l . the Feasible Real RAM model m or 
significance arithmetic [l8[. Implementations are for example MPFI ^sj, the 
iRRAM [21] and Mathematica 3l| respectively. 

All these methods have the same theoretical background. They are all 
practical instances of the model of Computable Analysis js^, 24, 16] used in 
computer science. While the Feasible Real RAM model directly implements 
the theory of Computable Analysis, the other mentioned methods all have 
their background in scientific computing. Looking closer at the various vali- 
dated methods in use, they all have in common implementing some kind of 
intervals for representing real numbers numerically. Therefore, the starting 
point here is looking at interval arithmetic for computing orbits For 
any time step n, let the phase point Xn together with its computational er- 
ror be represented by two fioating-point numbers and x" {x^ > x^) with 
given length m„ of the significand, called the precision, forming an interval 
[x^, x^]. The interval is an enclosure of the real value x„, that 
for all n. The interval length (i„ := xj^ — x^ gives a measure of the uncertainty 
about Xn and is therefore a kind of error. Interval arithmetic often models 
quantities which are not known exactly. But here, the true orbit (x„)„ can 
be, in principle, calculated to any given accuracy. Thus, in the present set- 
ting, the true object of interest is not an interval, but an approximation x„ 
of Xn together with an absolute error e^. The interval is only used for math- 
ematical convenience. To transform the interval to a fioating point value x„ 
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a;„:=rci(^^4^,mj (1) 



of precision m„, just do 



2 

where rd{x, m) performs a rounding to some floating-point number of pre- 
cision m nearest to x. Note that rounding to nearest is not unique if x 
is equidistant from two floating-point numbers. The absolute error e„ := 
\xn — Xn\ of Xn cau bc estimated via the interval length dn by 

en < ]^dn + r„ (2) 

where r„ is an error caused by the rounding operation rd{., .) in Equation 
(II]). An upper bound on r„ will be discussed later, for now it suffices to say 
that in general it is small compared to dn- 

The aim now is to calculate, for given initial value a; = Xq, G N and 
p G Z the orbit up to time N with relative error at most 10"^. That is, for 
(5'n)o<n<Af it should hold 

en = \Xn - Xn\ < IQ-^Xn < IQ-^ (3) 

Why using here and in the following the relative error and not the absolute 
error is discussed in some detail at the end of Subsection 13.21 The minimal 
m, fulfilling the precision requirement (|3]) on the relative error of Xn, which 
depends on x, N and p, is denoted by niminix, N,p). Now, a central quantity 
of this work is introduced, which is some complexity measure. Consider the 
growth rate of rriminix, N,p), 

I . y mmin{x,N,p) 
a[x,p) = hmsup — . 

N^oo N 

The loss of significance rate cr{x), which may depend on the initial value x 
is given by 

(7{x) = lim a{x,p). 

This quantity describes the limiting amount of significant precision being 
lost at each iteration step in the limit of infinite output precision. Significant 
means here the part of the digits being correct. A general treatment of 
this complexity measure is given in the next section. Roughly speaking, 
\a{xo,p)N + p-\d{10)] is the precision for any fioating-point number needed 
in an algorithm doing the iteration starting with xq and calculating to x^, if 
the output should be precise to at least p decimal places. Here, ld(.) is the 
logarithm to base 2. 
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2.1. Dynamic behavior of the logistic equation 

Before analyzing the different numerical behavior, it is worth having an 
analytical look at the dynamical behavior of the system. Despite the fact 
that these results are well known 13, 0], they are reviewed here for the sake 



of self containment. 

First have a look at the fixed points of the logistic equation and their 
stability. In the range D = [0,1], the equation possesses exactly one fixed 
point x° = if /i G (0, 1] and exactly two fixed points x° = and x^^^ = 1 — ^ 

if /i G (1,4]. Looking at the derivatives f'^{x°) = /i and f'^ix^'^^) = 2 — /i 
gives the stability of the fixed points. Since \f'fj,{x°)\ < 1 for G (0, 1) and 
> 1 for /i G (1, 4], x° is a stable fixed point, an attractor for fi G (0, 1) 
and an unstable fixed point, a repeller for yU G (1, 4]. If /i = 1, the only fixed 
point x° is hyperbohc, that is |/i(a;°)| = 1. At /i = 1, a bifurcation occurs. 
If yU G (1, 3), x° becomes unstable and the newly occurring fixed point x^'^^ is 
stable. At /i = 3 a second bifurcation occurs and for /i > 3 both fixed points 
are unstable. 

Second, examine the basin of attraction of the stable fixed point. If /i G 
(0,1), the contraction mapping principle directly gives limn^oo fji{x) = x° 
for all X G [0, 1]. If /i = 1, observe that fi{x) < x holds for all x G (0, 1). 
Hence, any sequence (/"(a;))^, x G (0, 1), is strictly decreasing and bounded 
from below. So, also lim„_>.oo /"(x) = x° holds for all x G [0,1]. Finally, 
in the case jj, G (1,3), lim„_>oo /^^(a;) = x'^^^ holds for all x G [0,1]. For a 
proof, the interested reader is referred to the literature: 0], Proposition 5.3 
in Section 1.5. 

Finally, for /i > 3 the system goes into a region showing periodic behavior 
with period doubling bifurcations. Finally, for some /i < 4, chaotic behavior 
is reached. 

This analysis shows that in the parameter range /i G (0,3), the orbit 
tends to the stable fixed point for any initial value x G [0, 1]. Furthermore, 
there exists some closed interval I ^ D, which depends on n, containing the 
stable fixed point such that f^{I) C / holds and is a contraction on /. 
Next have a look at the computational effort in the various control parameter 
ranges. 

2.2. Numerical analysis of the computational complexity 

The logistic equation is implemented in various forms using an arbitrary- 
precision interval library. For that purpose, the already mentioned interval 
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library MPFI based on the arbitrary-precision floating-point number library 
MPFR, both written in C, is used. For each control parameter /i ranging 
from 0.005 to 4 and a step size of 0.005, the orbit for initial value x = 0.22 
is calculated up to iV = 2000. For each fi, the minimum precision rrimin 
needed to guarantee e„ < 10~^a;„ for = 0, . . . , is searched. Then, aest '■= 
"^min/N is calculated. First, is implemented using a natural interval 
extension based on the expressions fix{l — x), — x^) and ^ — — |)^. 
The natural interval extension is obtained by replacing any occurrence of the 
variable x in the expression by an interval x |26|. The results are shown in 
Figures [H [2] and [3] respectively. Second, the logistic equation is implemented 



using a centered form, actually the mean value form 26|, ll7|: Ti^^(a;) = 
/p(mid(a;)) + f'^{x){x — mid(a;)) where x is an interval and mid(a;) is the 
midpoint of x. The result is shown in Figure HI In the following, these 4 
calculations are referred to as 1 to 4 respectively. 

The interval computation is in agreement with the dynamical picture only 
in Calculation 4. While for fi G (0, 1), the results shown in Calculations 1, 2 
and 4 are in agreement with the dynamical analysis, 3 is not since it would 
suggest an exponential divergence of initially nearby orbits which is not true 
in reality. A similar situation occurs for /i G (1,3). Here, the Calculations 
3 and 4 are in agreement with the dynamical picture, 1 and 2 on the other 
hand not. The picture does not change if > 3 and hence it can be said 
that in the range /i G (1,4], the Calculations 3 and 4 are in agreement with 
the dynamic picture, while 1 and 2 are not. How can this be explained? 

2.3. Investigating Calculation 1 

This subsection deals with the explanation of the curve obtained by Cal- 
culation 1. For doing an error analysis of the logistic equation analytically, 
some idealizing assumptions have to be made. Generally, executing the it- 
eration in interval arithmetic, two types of error are present. First, error 
propagation solely due to the iteration and second the newly added rounding 
error caused by the calculation of In the following, only the error prop- 
agation is regarded. This means that there is only one primary made error 
caused by rounding the initial value to some floating-point number 

of some specified precision m. The next idealization is that the value of ix 
is assumed to be given with such a high precision that no interval represen- 
tation is needed. Finally, the value of in Equation ([2]) is neglected. The 
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recursion relation then reads 

with the interval length dn given by the recursion relation 

dn+l /^(-^n ■^n'^n •^n ~^ •^n-^n) 

= jJ'dn 

with the obvious solution dn = fi"'do. The absolute error of Xn according 
to Equation ([1]) can be bounded from above by 

Cn < ^dn = ^/i"cio- (4) 

The ideal assumptions require the somewhat unreal setting that the precision 
has to be set to some finite, but big enough value m for representing xq 
and a virtually infinite value moo for doing the iteration. To get a relation 
connecting m and the output precision p in (E]), some upper bound on d^ is 
needed. The value of do is given as the rounding error by representing xq as 
a floating-point number of precision m. For that, the well known estimate 

do < 2-"'+'xo < 2-™+^ (5) 

exists. Combining ([3]), (j4]) and (|5]) gives as a sufficient condition 

for n = 0, . . . ,N. So, the sufficient condition gives an upper bound on 

N,p) by 

mrmn{x,N,p) < ■ Id(lO) + ■ max(0, ld(/i))l . 

This finally leads to an upper bound for the loss of significance rate, 

< cr(x) < max(0, ld(/i)). 

The curve in Figure [1] shows that (Test exceeds the estimated bound 
max(0, ld(/i)) only slightly. So, the above made ideal assumptions seem to 
be valid. In j2lj, the logistic equation was also investigated for jj, = 3.75 
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using the exact real arithmetic package iRRAM. In the paper, the maxi- 
mum bounding precision needed to guarantee the correctness of the first 6 
decimal places are reported up to = 100000. Relating this quantity to 
nimin shows fuU agreement with the simulation results performed here. So, 
for > 1, the interval length dn increases exponentially in time n which is 
in contrast to the dynamic behavior for /j, e (1,3). The reason is that the 
natural interval approach implicitly, due to the dependency problem, takes 
account only of the global behavior of in the form of a global Lipschitz 
constant max{|/^(x)| : x G D} = /i. However, a local Lipschitz constant 
max{|/^(a;)| : x G governs the real error propagation at time step n 

and also describes the dynamic behavior. 




Figure 1: Estimated loss of significance rate for the logistic equation, formula — x). 
The crosses indicate the theoretical curve from error analysis. 
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2.4- Investigating Calculation 2 

Calculation 2 is similar to Calculation 1. An analogous analytic approach 
as in Calculation 1 gives the recursion relation 

and hence 




Figure 2: Estimated loss of significance rate for the logistic equation, formula — x^). 
The crosses indicate the theoretical curves from error analysis. 



The recursion relation fulfills therefore dn+i > fJ^dn and dn+i < Sfidn- As 
a consequence the bounds /i"(io < c^n < (3/i)"(io are obtained. In analogy to 
Calculation 1, an upper bound for rriminix, N,p), 

ruminix, N,p) < \p ■ Id(lO) + N ■ max(0, ld(3/x))l 
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and hence 

(t{x,p) < (t{x) < max(0, ld(3/i)) 

is calculated. A brief look at Figure |2] shows that this upper bound is too 
rough. On the other hand, dn > fido suggests that the bound derived in 
Calculation 1 is a lower bound, hence max(0, ld(/i)) < a{x) < max(0, ld(3/i)). 
This is actually verified by numerical evidence. 

2.5. Investigating Calculation 3 

Calculation 3 is explained here in the parameter range /i G (0, 1), where 
it is not in agreement with the dynamic picture. Nevertheless it should be 
mentioned that the natural interval extension used here seems to be in full 
agreement with the dynamic picture in the parameter range G [1,4] as is 
suggested by Figure [31 The curve seems to be identical to Figure H] in the 
range /i G [1, 4]. 




0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 



Figure 3: Estimated loss of significance rate for the logistic equation, formula ^ — /i(a;— i)'^. 
The crosses indicate the theoretical curve from error analysis. 

To explain the observed behavior, first note that for /i < 1, //^(x) < jix 
follows for all x G [0,1]. Hence, x„ < /i"xo holds for all n G N and the 
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orbit {xn)n tends exponentially fast to zero. A brief look at the expression 
of ffj,{x), f^{xn) = 4 ~ At(a^n — \Y i shows that the value of the second term 
in the difference tends exponentially fast in n to the value of the first term. 
Hence, the big values of the loss of significance rate for small values of can 
be explained by cancellation. So, the behavior may be explained solely by a 
typical phenomenon of floating-point arithmetic and not an effect due to the 
dependency problem in interval arithmetic. To give an analytical description 
of the problem, it is easier now to chance from interval notation to classical 
error analysis notation. 

First note that even xq may differ from the initial value x since the con- 
version to a floating-point number may cause the very first rounding error. 
Next, already mentioned, in calculating the orbit two types of error are 

present. First, error propagation due to the iteration scheme and second the 
rounding error caused by the calculation of Now, let Xn for some n G N 
be given. Then the true error after one iteration step is — Xn+i- Since 
in reality not f^{xn) is calculated but some erroneous approximation 
the true error can be written as Xn+i — Xn+i = ffi{xn) — f^i{xn)- Inserting a 
constructive zero gives a sum 

Xn+l - = ift^iXn) " //.(a^n)) + (//.(^^n) " //.(^^n)) (6) 

of two terms. The first term describes solely the error propagation while the 
second term gives exactly the newly produced error due to the approximate 
calculation of 

Let us fix some n G N and consider the absolute error of Xn+i- To get the 
formulas more compact, set g^{x) := fi{x — |)^. Then, 

\Xn+l - Xn+l\ = I ((f) - g,,iXn)) " (f " ^^(^n))! 

< l((?)^(^n)) - ((?) - ^m(x„))| + 1(1) - f I 



4; ^^i\-^njj \\4y ^fi\-^njj\ ~l" 1(4) 4 I 



follows. Let m be assumed to be the actual precision under calculation at 
time n. The last term in the previous inequality can be estimated the fol- 



lowing way. As discussed in [35| , the rounding error produced in calculating 
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g^j, can be estimated by 

\9^.ixn) - 9^ixn)\ < lMK2-"'\g^{xrr)\ (7) 

where K is the number of rounding operations performed in computing g^. In 
the case considered here, K = 4 follows. It is further crucial to mention that 
the factor 1.06 is only valid if < 0.1 ■ 2"^ holds so that the precision must 
not be chosen too small. Furthermore, with \g^{xn) — gfi{xn)\ < — Xn\ 
it follows 

\xn+i - Xn+i\ < 2-"^|(f) - g^{xn)\ +2-^-11 + lMK2-^\g,{xr.)\ 

-\- jJj I X^i I 

< 2-™(|(f)| + \g^{xn)\ + f + imK\g^{xn)\) + - Xn\ 

< 2— ((1 + 2—)^ + (1 + 1.06if2-™)|(7^(x„)| + ^ 
+ lMK\g^{xn)\) + /u|x„ - Xn\ 

< 2-™f (1 + 2"™ + 1 + 1.067^2-™ + 1 + 1.06ir) 

~l~ A'' I '^n '^n I 

< C/i2~™ + - Xn\ 

where C > holds. In other words, one obtains the recursion relation 
Cn+i < /ie„ + Cyu2~™'. Iterating the recursion gives e„+i < C'/i2~'" ^^^q /x'^ + 
u"+ieo < C-^2-™ + u"+i2-'"xo. 

As already mentioned, x„ is bounded from above by x„ < /^^xq. To come 
to a sufficient condition for the precision, also a lower bound is needed. First 
observe that /^(x) > yux(l — a) holds for all x < a, a, x G [0, 1]. Hence, for 
/i < 1, x„+fc > /i"'Xfc(l — Xfc)" follows. This gives the sufficient condition 
(j_E_2-m ^ ^n+i2-m^^ < 10-P/i"^i-'=Xfc(l -Xfc)"+^-^ < 10-Px„+i 

on the precision. Note that n + 1 > A; > 0. Then, an upper bound on 
m^i„(xo, N,p) is given by 

m™„(x, N,p) < \p ■ Id(lO) + {N- A;)(ld(i) - ld(l - x^)) + C] 

with C = Id^Cj^ + /i^Xo) — Id(xfc). This leads to an upper bound on the 
loss of significance rate given by a{x,p) < Id(^) — ld(l — Xk) for all G N. 
Since x^ — )■ follows for — )• oo, the final result on the loss of significance 
rate is 

a(x,p) <a(x) <ld(i). 

The curve in Figure [2] shows that this upper bound is in full agreement with 
the numeric result. 
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2. 6. Investigating Calculation 4 

The observation at the end of the subsection describing Calculation 1 
directly leads to the already introduced mean value form. The calculation is 
shown in Figure |H This calculation is the optimum of both, Calculation 1 
and 3. The curve reflects in the parameter range /i G (0, 3) well the dynamic 
behavior. 



1.0 
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0.4 
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0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 

Figure 4: Estimated loss of significance rate for the logistic equation, meanvalue form. 

Furthermore, in the range /i G [3, 4], the curve suggests a relation between 
the loss of significance rate and the Lyapunov exponent A(x) for the logistic 
map: 

^(a;) = j^max(0,A(x)) 

for all n G (0,4]. For a curve of the Lyapunov exponent of the logistic map 
see jsl . This relation will be shown in the next section for general dynamical 
systems on the interval. Furthermore, it will be shown that the algorithm 
based on Calculation 4 is optimal in some sense. 



13 




But before, some crucial reflections governing the analysis in the next 
section. The mean value form representation, on which the calculation is 
based, can also be seen from a different viewpoint. Have again a look at 
Equation ([6]). The true error is the sum of the error propagation (first term) 
according to the iteration and the rounding error due to the computation 
of (second term). The first term of Equation can be handled using 
the mean value theorem, |/;,(x„) - /^(x„)| = |//,(?/n)| ■ \x„ - x„| with ?/„ e 
[xn — Cn, Xn + e„] . This gives directly the bound 

\ff,{.Xn) - ft,{,Xn)\ < SUp(|/;^([x„ - e„, + en])|)e„. 

The second term can be estimated in a similar way as was done in ([7]) by 
\f^.{in) - /M(*n)| < lMK2-^\f^{x^)\ 

where K = A because there are 3 arithmetic operations and the rounding of 
ji. Using the fact that f^{x) < ^ holds and /^(a;) < a; if /i < 1, the unknown 
value \f^{xn)\ can be estimated from above. This calculation shows that 
there exists a recursive equation on an upper bound e„ on e„ for all n: 

en+i = L(xn,e„)e„ + 1.06ir2-™E^(a;„), eo = 2^"^ (8) 

with L(x, e) := sup(|/^([a; — e,x + e])|) and 



E^{x) :-- 



X if /i < 1 
if /i > 1 
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This description, which is in line with the analysis of Calculation 3, is equiv- 
alent to the interval description using the mean value form. Instead of using 
intervals, pairs of the form value x„ and corresponding guaranteed error 
bound e„ is used. This approach is an automated error analysis called run- 



ning error analysis From a technical point of view, the representation as 
value and error has the advantage that the rounded values x„ are calculated 
as usual in floating-point arithmetic except that arbitrary-precision floats are 
used. The guaranteed error bounds may be calculated using interval arith- 
metic according to ([8]), to really guarantee a validated bound. Only a fixed 
precision is needed for calculating the error bounds. Similar results as in 
Figure H] are reported in [2| by using a method analog to the one presented 
here [sj. However, the connection to the Lyapunov exponent is not made in 
i. 
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Before continuing, three remarks. First, interval libraries are primarily di- 
vided int two types concerning their representation of an interval 29|| : There 



exist libraries using the infimum-supremum representation of intervals, like 
MPFI, and there exist libraries using the midpoint-radius representation of 
intervals. If arbitrary precision is needed, the inf-sup libraries have the dis- 
advantage that two floating-point variables with high precision are needed to 
represent an interval. Contrary to that, like the value and error description, 
in mid-rad libraries only the midpoint of the interval needs a high precision 
floating-point variable. The radius can be stored in a floating-point variable 
which need not have a high precision. Clearly, the mid-rad concept has an 
computational advantage in the case considered here over the inf-sup con- 
cept. But the dependency problem of interval arithmetic persists. Second, 
also the iRRAM package implements mid-rad intervals and has therefore to 
cope with the dependency problem. However, it also permits an optimized 
way for computing the iteration based on a similar algorithm as described 



above [22|. Third it should be mentioned that, executing the first three pre- 
sented calculations in Mathemathica using significance arithmetic, exactly 
the same results are obtained. This shows that also significance arithmetic 
suffers from the dependency problem as interval arithmetic does. This is 



already noted in [32 



3. The general algorithm and its complexity 

Let D be a compact real interval and f : D ^ D a. self mapping. In 
the following, / is assumed to be continuous on D, two times continuously 
different iable on D and /" is bounded. Furthermore, / and /' are assumed 
to be computable in the sense of Computable Analysis. The definition of a 
computable real function is given below. 

In this section, a general algorithm for computing the iteration 

Xn+l = f{Xn), XqG D (9) 

is presented. To be more precise, for given x & Q, N E N and j9 G Z, this 
algorithm computes a finite part {xn)o<n<N of length of the true orbit 
(x„)„gN with initial value xq = x. Each computed value x„ of this finite 
trajectory has a relative error of at most 10~^: | I < 10-P\xn\ for 

all n = 0,1, . . . , N. The correctness of the algorithm and its relation to 
Computable Analysis is shown. Finally, its complexity is examined. 



15 



3.1. Computability issues and specifying the algorithm 

The set of all computationally accessible real numbers are the floating- 
point numbers of arbitrary precision and arbitrary exponent range denoted 
by M. A floating-point number is a real number of the form x = s- 2^~* where 
t G N is the precision, e E Z the scale and s e Z where |s| G {0, 1, . . . , 2* — 1} 
is called the significand. To get a unique representation of x for given t, 
\s\ > 2*~^ is assumed if x 7^ and e = if x = 0. Since actually no bound is 
assumed on the precision and the scale, the set M C M is the set of the dyadic 
real numbers and therefore countable infinite. Thus, M forms a natural basis 
for computability considerations over finite objects. Consider some floating- 
point number x G M, then the scale and the precision are two properties 
of different type. While the scale is a direct function of the value of x, the 
precision is clearly not. Reversely, let x G M be some real number and x G M 
a floating-point number representing x. Then the scale of x is generally 
determined by x while the precision can be chosen arbitrary. Regarding x as 
a data structure, then x has as its essential property the precision. In object 
oriented notation, the precision of x can be written as x.t. 

Any real number x is represented in an algorithm concerning numerical 
computation by a pair x consisting of a floating point number x.fl G M 
of arbitrary precision x.fl.t approximating x and a floating-point number 
x.err G M of fixed precision giving an upper bound on the absolute error, 
\x.fl — x| < x.err. Reversely, any such pair a? G can be seen as the real 
interval [x. fl— x.err, x.fl+x. err]. If x G [x.fl — x.err, x.fl+x.err] holds for 
some X G M, then x is called an approximation of x. To represent a single real 
number, a sequence {xn)n£n of such pairs x are needed. A sequence {xnjnevi 
is called a floating-point name of a real number x, if any Xn approximates 
X, lim„^oo Xn.fl = X, lim„^oo x^.fl.t = 00 and lim„^oo x^.err = holds. 
Clearly any real number has a floating-point name. 

As already indicated, it is a straightforward task to define what a com- 
putable function / : M — )■ M is by using classical computability theory 
over finite objects. Additionally, computability over integers, computabil- 
ity of functions with mixed arguments and computable predicates are de- 
fined in the same manner [33|. Consider a function f : D ^ D, D <Z M. 
and a pair / of two functions f.fl : M — )■ M and f.erf : — )■ M hav- 
ing the following property. For any approximation x of some real num- 
ber X G -D, the pair f{x) = {f.fl{x.fl),f.erf{x)) is an approximation of 
/(x). Thus, f.erf gives an upper bound on the absolute error of f .fl{x.fl), 
\f.fl{x.fl) — /(x)| < f.erfix). Considering / as an interval function, the 
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above property is just the fundamental property of interval arithmetic, [27 
Property 2.12. Then, / is called an approximation function for /. Now 
consider an approximation function / for / such that for all x G -D and 
any floating point name {xn)nm of 2;, {f{xn))nen is a floating-point name 
of f{x). Such an approximation function is called approximation- continuous. 
Additionally, if the two functions /.// and f .erf of an approximation func- 
tion / are computable, then / is called a computable approximation function. 
Finally, f : D ^ D is called computable, if there exists a computable approx- 
imation function / for / which is approximation-continuous. 

The algorithm with the specification described at the beginning of this 
section reads 

1 Input parameter: x, N, p 

2 Initialize precision m 1 

3 do 

4 Initialize value and error x ■(^ rd{x,m) 

5 f orn = to iV do 

6 If prec{x,p) = true then 

7 If not printed print n, x.fl, x.err 

8 else break 

9 x^ fix) 

10 end for 

11 m 4- m + 1 

12 while prec{x,p) = false 

where / is an approximation-continuous approximation function for / spec- 
ified below. To initialize x, a rounding function re? : Q x N —)■ is needed 
where rd{x, m).fl is a floating-point number of precision m being the exactly 
rounded value of x for some rounding convention, in the following nearest. 
Clearly, the value rd{x,m).err is an upper bound on the absolute rounding 
error, rd{x,m).err = ^ulp{x) if the rounding mode is nearest. The predicate 
prec : X Z — i- {true, false} is a test whether the relative error of x.fl is 
bounded by 10^^. The semantics reads: 

If cc G approximates x G M and prec{x,p) = true holds, ^^^^ 
then \x.fl — x\ < 10~''|x| follows. 

While the object oriented notation is convenient for a compact and in- 
structive description of the algorithm, in the following analytical analysis an 
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abbreviation for this notation is sometimes more handsome. As in the hne 
of the preceding section, floating-point numbers and functions are indicated 
by a hat: x := x.fl and / := /.//. An over-bar indicates an error bound: 
e := x.err and erf := f.erf. Hence, x is equivalent to (x,e) and / is 
equivalent to {f,erf). 

Finally a remark on optimization. The algorithm is not optimized in 
performance. Including performance issues, in Line 11 something like m 
a ■ m + b can be used where a > 1 and c e N are constants. Here, the aim is 
to find the minimal m to guarantee some given upper bound on the relative 
error of x„. 

3.2. Computability and correctness 

It is clear that the rounding function rd is computable. So let us begin 
with the predicate prec. 

Proposition 3.1. The predicate 

prec(x,p) := | " TTIo^l^-^^l (H) 

[ false otherwise 

is computable and satisfies / flOj) . 

Proof. Let x be an approximation of x. If x.err < -^p^^\x.fl\ holds, then 
x.err < 10~P(|a;.//| — x.err) follows. Using \x.fl\ < \x.fl — a;| + \x\ < 
x.err + \x.fl — x| < x.err < 10~^(|a3.//| — x.err) < 10^^|x| follows. 

The predicate (fTTl) only uses the approximation x, basic arithmetic and 
finite tests. Hence, this formula is computable. □ 

Note that the definition of the predicate also gives true in the singular 
case where x.fl = and x.err = and hence x = 0. 

An algorithm for computing /.// is possible by assumption. To derive 
an algorithm for computing f.erf on the absolute error, return to Equations 
dSD and dH]). 

Proposition 3.2. Let x E D he given and x an approximation of x with 
x.fl G D. Assume that f.fl{x.fl) computes the value f{x.fl) up to a cor- 
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rectly rounded last bit in the significanda Furthermore assume f.fl{x.fl).t = 
x.fl.t. Then the absolute error of f{x) is bounded from above by 



L{x) ■ x.err + 2-^-^'-* ■ \f.fl{x.fl)\. (12) 
Here, L(x) = sup(|/'([a;.// — x.err, x.fl + x.err] n D)\). 

Proof. Equation ^ gives \ f{x) — f{x)\ < \ f{x) — f{x)\ + \f{x) — f{x)\. Using 
the mean value theorem, — < sup(|/'([a; — e, x + e])|)e follows. 

According to the assumption on / and Theorem 2.3 of [l^], \f{x) — f{x)\ < 
2~™'|/(ai)| holds where m is the precision of x. □ 

Corollary 3.1. Let f be as specified in the beginning of this section, f.fl 
and L{x) specified as in Proposition \3.S[ Then there exists a function L{x) 
with L{x) < L{x) < Lmax for some L^ax > such that f with f.erf{x) = 
L{x)-x.err+2^'^-^'--^-\f.fl{x.fl)\ is an approximation- continuous, computable 
approximation function of f . 

Proof. Let L{x) be some computable upper bound of L{x). L{x) can be 
computed by global optimization, for example by using interval arithmetic. 
Since /' is continuous and D compact, L{x) is bounded. So, L{x) < Lmax 
for some L^ax > 0. Also, / is computable. Using Proposition 13. 2[ it follows 
that / is also an approximation function of /. Remains to show that / is 
approximation-continuous. Let {Xn)n be some floating-point name of x G -D. 
Clearly limn^oo f -fK^n-f I) -t = oo holds. Since lim„_>.oo err = and the 
sequences {L{xn))n and {\f.fl{xn.fl)\)n are bounded, lim„_^oo f.erf{xn) = 
follows. Furthermore, by this result and the statement of Proposition 13. 2[ 
also lim„^oo f-fK^n) = x holds. □ 

To summarize, the iteration (jO]) is performed in the algorithm by iterating 
a value x„ approximating Xn with an upper bound on its absolute error e„ 



^This assumption is pragmatic. The already mentioned software package MPFR 
implements this specification. The problem of achieving this task for transcenden- 
tal functions may be of unknown cost and is known as The Table Maker's Dilemma, 
see http : //per so . ens- lyon. fr/ jean-mi chel .muller/Intro-to-TMD .htm. Addition- 
ally note that this assumption can be weakened without abandoning the main statements 
of this work. 
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according to 



Xo = rd{x, m) (13) 
eo = 2-™|xo| (14) 

where I/(x„, e„) is a computable upper bound on L(x„, e„) as described in the 
preceding corollary and m the precision of any floating-point number involved 
at that stage. This is Line 9 in the inner f or-loop of the algorithm which 
is executed with successively increasing precision m, controlled by the outer 
do-while-loop. Finally, it has to be shown that this outer loop eventually 
terminates. 

Proposition 3.3. Let x & D with x ^ be given and {x„i)m>i a floating- 
point name of x obeying Xm-fl-t = m. Then \im.m^ooP^Gc{xm-,p) = true 
follows for all p G Z. 

Proof. Since x ^ and limm_5.oo Xm-err = 0, there exists some M G N such 
that for all m > M, < \Xm.fl\ and x^.err < 2{i+iq-p) \^\ holds for all 
m > M. Then prec{Xm,p) = true for all m > M. □ 

The next proposition makes the link to Line 9 in the algorithm. 

Proposition 3.4. Let Xn be the n-th element of the orbit of the recursion ^ 
and {{Xn)m)m>i 0- sequence given according to the recursion equations [W^] 
and ( (771 ) '"^^^^ increasing precision {xn)m-fl-t = m. Then {{xn)m)m>i is a 
floating-point name of Xn- 

Proof. Let Lmax according to Corollary 13. II and M > sup{|x| : x E D} such 
that \xn\ < M holds for all n. Then Equation f|T^ leads to e„+i < Lmax^n + 

2-M. Iteration gives e„ < lL.eo+2-M E^lJ ^l. < S-^MELo^L.- 
Hence, for n fixed, limm__>oo(^c„)m-err = follows and consequently also 

iim^_>o53 (£Cjj)n2.y/ Xfi- n 

These two propositions finish the correctness proof of the algorithm. They 
show that, if Xn ^ for n = 0, . . . ,N, the outer loop eventually terminates 
for any p G Z. 

The drawback of the algorithm is, that in the = for some 

n < N, the computation does not terminate. This is only due to the fact 
that the relative error controls the outer do-while-loop. If the absolute error 
would be used instead, this drawback is eliminated. However, controlling the 



Xn+l — fi^n) 
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relative error is more general. Consider for example a dynamics with positive 
phase space, the concentration of a substance for example. If the value varies 
in time over a wide range in scale, it is fortunate to illustrate the orbit in 
a logarithmic plot. If the relative error is controlled, the error bars in the 
plot are constant, in contrast to large varying error bars in the case where 
absolute errors are used. 

Absolute errors arc in the line with Computable Analysis. Replacing the 
test prec{x,p) by the test on x.err < 10^^ in the algorithm would give a 
segment (a;„)o<n<Af of the orbit with accuracy \xn — < 10^^. It is now 
straightforward to see that the function g : D xN ^ D with g{x, n) := f'^{x) 
is computable. Here, a function : D x N — )■ D is computable if there 
exists a computable approximation function : x N — )■ for g which is 
approximation-continuous with respect to the first argument. 

3.3. Computational complexity 

After having presented the preliminary work, the main issue of the paper 
is addressed - the computational complexity of the presented algorithm. The 
complexity measure of interest here is the loss of significance rate already 
introduced informally in the previous section. Here is the formal definition. 

Definition 3.1. The minimal precision, for which the described algorithm 
eventually halts is denoted by m^j„(a;, N,p), where x, N and p are the cor- 
responding input parameters. The growth rate of mmin is given by 

a[x, p) hm sup . (15) 

Then, the loss of significance rate cr : Q fl D ^ M is defined by 

(t{x) := lim a{x,p). (16) 

p— ^-oo 

To achieve bounds on the loss of significance rate, the drawback of the 
preceding subsection also makes problems here. If = for some n G N, the 
loss of significance rate may be unbounded. Therefore, one more assumption 
in addition to the ones on the dynamical system stated in the beginning of 
this section has to be made. 

Assumption 3.1. The dynamical system (-D, /) is assumed to have the 
properties already mentioned in the beginning of this section and further- 
more, for any orbit {xn)n under consideration, x„ 7^ holds for any n e N 
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as well as 

j.^ ld(min{|x„| :n = 0,l,...,N}) ^ ^ 

N^oo N 

If only a finite range in scale is relevant, the additional assumption is no 
loss of generality. An example is the logistic equation where G -D but has 
no distinguished role. Instead of considering {D,f), consider the following 
dynamical system {D, /). Choose some M > — mm{D) and set D := {x + 
M \ X e D} as well as f{x) := f{x - M) + M for all x E D. Then (D, /) 
fulfills the additional assumption. Furthermore f'{x) = f'{x — M) holds and 
therefore there is no substantial difference in the complexity analysis of the 
algorithm between the original system and the modified system. 

First, the boundedness of a{x) is shown. 

Proposition 3.5. Let {D,f) be as in Assumption \3. 1\ and m^,„ (x. N. p) as 
in Definition \3.1[ Then, for given p G Z, there exists a constant C > 0, 
depending on f , such that mjnin{x, N,p) < C ■ N + o{N) holds for all N eN, 

X eQnD. 

Proof. According to the requirements made on {D,f), there are some con- 
stants L > I and M > such that e„+i < Le„ + 2"'^M holds for all 
n G N and all precisions m. Analogous to the treatment in the proof of 
Proposition [331 iteration gives e„ < 2"™MX;Lo^*' = 2-"^M ^7-i~^ - Let 
B{N) := min{|x„| : n = 0, 1, . . . , A^} > 0. Then, for all n < N, e„/|x„| < 
en/B{N) < (^3^2-L^+i follows. If now (xr3^2-L^+i < ^ 

holds, prec{{xn, e„),p) = true for all n = 0, . . . , A^. This leads to the bound 
mminix, N, p) < ld(L) ■ A^ + max(l, Id(f^) - ld{B{N)) + p ■ Id(lO) + ld(l + 
10-P)). □ 

Corollary 3.2. Let {D,f) be as in Assumption \3.1[ (j{x,p) as in fTE\) and 

a{x) the loss of significance rate. Then, for given p G Z, there exists some 
constant C > such that a{x,p) < a{x) < C holds for all x E Q (1 D . 

In the following, the main statements of this paper are be formulated: A 
lower and an upper bound for the loss of significance rate is given. Further- 
more, the relation of these bounds to the Lyapunov exponent X{x) is shown. 
Before the theorem is stated, for sake of completeness, the definition of the 
Lyapunov exponent and its basic properties are presented. 
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Definition 3.2. Let {D, /) be a dynamical system, D C R compact and 
f : D ^ D continuously differentiable on D. Then the Lyapunov exponent 
at a; G -D is defined by 

^ n— 1 

X{x):= lim -J2H\f{f{m (17) 

71,— ^no n ' ' 



k=0 



if the limit exists. 



The Lyapunov exponent may depend on x. However, the following prop- 
erties hold: 

(a) If {D, f) has an invariant measure p, then the limit in Equation (1T7|) 
exists p-almost everywhere. 

(b) Furthermore, if p is ergodic then A(x) is p-almost everywhere constant 
and equal to 

' \n{\f'{x)\)p{dx). 



'D 

These properties are a direct consequence of the Birkhoff ergodic theorem, 
see 



15|, Theorem 4.1.2 and Corollary 4.1.9. Now the first theorem. 



Theorem 3.1. Let {D,f) be as in Assumption \3.1\ cr{x,p) as in l[TS\) and 
X{x) the Lyapunov exponent of {D,f). Then cr(x,p) > max(0, A(x))/ln(2) 
holds for all X E Q n D , p E Z, if A(x) exists. 

Proof. Let G N be given and M > a constant with < M for all n G 
N. According to Equation f lT^ and Proposition 13. 2 ^ e„+i > |/'(x„)|en holds. 

Iteration gives > |xo|2-™nl'o If'Ml So, |^ > ^ nl'o \fM\ 
follows. A necessary condition for the algorithm to terminate is therefore 
§^2— nl"o'l/'(a;n)| < This gives the bound on m™„(x,Ar,p) > 

E!L~o^ld(|/'(a;fc)|) + ld(^)+p-ld(10)+ld(l + 10-f). Following the definitions 
of <j{x,p) and the Lyapunov exponent, cr(a;,p) > A(x)/ln(2) follows. □ 

Before a realistic upper bound on <7{x,p) can be presented, one more 
definition is needed. 



Definition 3.3. Let a > then define a function rja : (0, oo) — M by 

ria{x) := 
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ln(x) if X > a 
ln(a;) if x < a 



Furthermore, for any a > define 

71—1 



Proposition 3.6. For all a > there exists some constant C > such that 
^a{x) < C holds for all x & D. Furthermore, if the Lyapunov exponent X{x) 
exists, X{x) < Xa{x) holds. 

Proof. According to tlie requirements made on {D, /), / is Lipschitz with a 
Lipschitz constant L > 0. Furthermore, let a > be given. Then for all 
n eN, ^T.k=lVa{\f{fHx))\) < ln(max(a,L)) holds. Hence it follows the 
upper bound on limsup^^^^Y.k=oVai\f'if^ix))\) < ln(max(a, L)). The 
second assertion follows from the fact that ln(x) < r]a{x) holds for all x > 0, 
a > 0. □ 

Proposition 3.7. Let x E D he given. If X{x) exists, then also the limit 

lim A„(x) =: A(a;) (18) 

exists and X{x) > X{x). 

Proof. Since ln(a;) < ria{x) < rj^^x) holds for all a; > 0, < a < /3, also 
X{x) < Xa{x) < A/3(x) follows. Letting a — )■ 0, a > 0, the assertion follows. 

□ 

Theorem 3.2. Let {D,f) be as in Assumption \3.1[ cr{x,p) as in ( 173]) and 

X{x) as in fT^) . Let x G Q fl D be given, then for any e > there is some 
Po G ^ such that for all p > po, 

c^ix,p) < j^max(0, A(a;)) +e 

holds if X{x) exists. 

So there is the following bound on the loss of significance rate. 

Corollary 3.3. Let {D,f) be as in Assumption \3.1\ a{x,p) as in l[T5\) . a{x) 
the loss of significance rate, X{x) as in l[T^) and X{x) the Lyapunov exponent. 
Then, 

j^max(0, A(x)) < o-(x) < max(0, A(a;)) 
holds for all x E Q n D if X{x) exists. 
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Before the proof of the theorem can be presented, the following lemma is 
needed. 

Lemma 3.1. Let e > and a > y/e. Then for all x > 0, 

ln(x + e) < r]a{x) + y/e 

holds. 

Proof. There is nothing to prove in the case e = 0. So let e > 0. Two cases 
are considered. 

1st case: x > a. Then the inequality reads ln(x + e) < ln(x) + y/e which 
is equivalent to x > — r^^=r- r- Since — r^^r-r < 4= < « < the assertion 
follows. 

2nd case: x < a. Then the inequality reads ln(x + e) < ln(a) + y/e 
which is equivalent to x < aexp(A/£) —e. A sufficient condition to prove the 
assertion is a < aexp(A/e) — e which is equivalent to a > ^^^f^^-^_i ■ This 
was already proven in the first case. □ 

Now everything is prepared to prove Theorem 13.21 

Proof of TheoremlKM Let N eN, B{N) := min{|x„| : n = 0, 1, . . . , A^} > 
and M > a constant with |a;„| < M for all n E N. Starting with Equation 



Mj) and iterating gives 



Define 



and 



n— 1 n n—1 

eoY[L{xi,ei) + 2-"'^\xk\ Y[L{xi,ei) 

1=0 k=l l=k 

n n—1 n n—1 

z 1*'^ I n ^ ^^2-"^ Yi n 

k=0 l=k k=0 l=k 

n n—1 

spn ■■= YYlL{xi,ei) 

k=0 l=k 



SP{N) := max{spn : n = 0,1, . . . , N}. 

Then, e„ < M2~'^SP{N) follows for all n < N. A sufficient condition for 
the algorithm to terminate is given by -^^2~^SP{N) < Hence, 

m^^x, N, p) < max(l, C - \d{B{N)) + \d{SP{N))) 
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follows with C = ld(M) + p ■ Id(lO) + ld(l + IQ-p). Using the Assumption 
13.11 leads to 

cr{x,p) < -^max (o,limsup^ln(S'P(A^)) ) . 
By definition, lim sup^^^^ iri(SPiN)) ^ sup^^^^ follows and hence 

(j{x,p) < 7-^—- max ( 0, lim sup ^^^'^ 



ln(2) V „->oo n 
Next let 

rt-l 

Pk,n ■= Y\_H^hei) 

l=k 

for G N and k < n, and furthermore 

P„ := max{pfc,„ : k = 0,1, . . . ,n}, 
then sj9,i < {n + l)Pn follows for all n G N. This gives 

In(Pn) 



(j{x,p) < max I 0, lim sup 



ln(2) V n-^oo n 
Let i^(n) G {0, . . . ,n} be the smallest number such that nr=r7^(n) L{xi,ei) = 
Pn- Then consider ^^^^^ = ;^ ^"J^j-^-, ln(L(a;;, e;)). Let V be a Lipschitz 
constant of /', then L{xn,en) < \ f'{xn) \ + L'2en holds for all n G N. Conse- 
quently, there exists some L > L' such that L{xn,en) < |/'(a^n)|+2Iv e„ holds 
for all n G N. This inequality leads to L{xn,en) < \ f'{xn) \ + 2L M ^|.°Q^p < 
+ 2I'm ■ IQ-P. Inserting gives 

i^<- V ln(|/'(x,)|+2l'M-10~^). 
n n ^-^ 

l=K{n) 

Now let e > and < a < 1 be given. Then choose po & N such 

that M ■ 10 < min(a, ln(2)|) holds. Then for all p > po, the above 
lemma gives 

E (r,.(|/'(xO|) + ln(2)|) (19) 

i=i<'(n) 

^ 71—1 

^M2)| + - 5^ (20) 

Z=_ft:(n) 
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Consider the sequence {K{n))n(^N- Observe that, first the sequence {K{n))nizfq 
is increasing and second if K{n + 1) > K{n) for some n G N, then K{n + 1) = 
n or K{n + 1) = n + 1. There are two cases. 

1st case: {K{n))nm is bounded. Then, there exists some constant Nq 
such that K{n) = K{No) holds for all n > Nq. Choose now a small enough 
such that Aq,(x) < \{x) + ln(2)| holds. Then, compute the upper limit 
to limsup„^^ i Er=K{n) Va{\f'{xi)\) = limsup„_^^ ^ Er=i^(7Vo) Va{\f'{xi)\) = 
limsup„^o^ i Yl'iZo Va{\f'{xi)\) = Xa{x). By taking the upper limit of ([20]), 
limsup^^^ < ln(2)| + Ac,(a;) < \n{2)e + \{x) follows. 

2nd case: {K{n))n£n is not bounded. Then, for any 6 > and any 
A^o ^ there is some n > Nq with ^"'•^"^ < 5. Since, by definition, 
J2]=o ^^{L{xi,ei)) < ln(P„) holds as well as \f'{xi)\ < L{xi,ei), the inequality 
^T.i=oH\f'ixi)\) < follows. This shows X{x) < 0. 

Next it is stated that for all e > and p > po, limsup^^go ^^^^^ — 1^(2)^ 
holds. This shows (t(x,p) <e = max(0, A(x))+£: < max(0, A(x))+£. 

Assume otherwise. Then, for some e > and N eN, first ^-^^^^ > ln(2)e 
holds and second A(x)-ln(2)| < ^ ln(|/'(x/)|) < A(x)+ln(2)| holds for 
all n > K{N). Using (120|) . the first expression gets J2h^K(N) Va{\f'{xi)\) > 
ln(2)|. Choose a small enough such that ?7a(|/'(x;)|) = ln(|/'(x;)|) holds 
for all / < N. Then, for sufficiently high p, jf^Z~KiN)H\f'{xi)\) > ln(2)| 
follows. In the second statement, the sum can be split the following way: 
jjEfJo^~'H\f'ixi)\) + jjElK\N)H\nxi)\) < A(x) +ln(2)|. The first 

addend on the left side is bounded form below by -^^(A(a;) — ln(2)|) > 
A(a;) — ln(2)|, the second addend is bounded from below by ln(2)|. Hence, 
A(a;) + ln(2)| < j^EZ'o^ H\f'ixi)\) < X{x) + ln(2)| follows, but this is a 
contradiction. □ 

In the end, it is shown that, if A(x) = A(x) holds, the algorithm presented 
here is optimal with respect to the loss of significance rate. This means that 
no algorithm with the specification presented at the beginning of this section 
has a lower loss of significance rate than the algorithm presented in this 
section. 

Proposition 3.8. Let {xn)n be an orbit of the dynamical system {D, f) and 
X{xo) > 0. Then, for any e > there exists an Nq E N such that for any 
N > No there is some 6 > such that the following holds. Let an initial 
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value uq G [xq — 5,Xq + S] H D be given and consider the corresponding orbit 
{yn)n- Then, 

holds. 

Proof. Let > be given. Then there exists some Nq e N such that 
jf^n=iM\fM\) > A(xo) - e' holds for all N > N^. For given iV > A^q 
there is some S' > with min{|/'(a;„)| : < n < A^} > S', otherwise A(a;o) 
would not exist. Consider now an orbit (y„)n with \yn — Xn\ < '^'^jj fo^ 
n — 0, . . . , N where L' is a Lipschitz constant of /' and < cr < 1 arbitrary. 
Then, for n = 1, . . . , A?", the following estimation holds. 

IVn - Xn\ = IfiVn-l) - f{Xn-l)\ 

= \f'iXn-l)iy„-l - + ^f"{Cn-l){yn-l " 

> (|/'(a;„_i)| - ^|/"(Cn-l)| • IVn-l - Xn-l\)\yn-l " Xn-l\ 

> {\f'{Xn-l)\ - ^L'ly^^i - Xn-l\)\yn~l " Xn-l\ 

> (|/'(a;n_i)| - a6')\yn-i - Xn-i\ 

where {n-i e [xn-i,yn-i] U Xn-i]. Iterating finally gives \yN - xn\ > 
Y[n=o(\f'(^n)\ — o'S')\yo — xo\- Now determine some constant C > such 
that ln(|/'(a;„)| - ad') > ln(|/'(a;„) |) - C holds for aU n = 0, . . . , iV the 
following way. A short calculation shows that this is equivalent to C > 
^^(l + M^F^)- Using ln(l+ |^,(4;_^,, ) < ln(l + ^) = ln(l + ^) < ^ 
finally gives C > as a sufficient condition. Set C — Let £ > be 
given. Set C = e' = e/2, then 

7V-1 ^„ iV-l 

EMI/'(^n)| - j-^) > J]ln(|/'(x„)|) - Arc 

n=0 n=0 

> N{\{xo) -e' -C) 

and hence 

It/AT — Xjvl^e \yQ — XQ\ 

follows. □ 
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Proposition 3.9. Let {D, f) be as in Assu'mption \3.1\ and x G Q fl D given 
such that X{x) exists and X{x) > 0. Consider an algorithm computing an 
initial segment (a;„)o<n<Af of the orbit of x with relative error < 10'^ for 
some e N, p e Z. Then the algorithm has a loss of significance rate 
crix) > 42)A(x). 

Proof. Let < e < X{x) be given and big enougli sucli tliat tlie previous 
proposition holds for some 6 > 0. Clioose some y E [x — 6, x + 6] H D , y ^ x . 
Let po be big enough such that 

2M10-P° < e^(^(^)'")|?/ - x\ (21) 

holds, where M > such that |x| < M for all x G D. Consider y as the 
initial value of another orbit {yn)n- Then, with the above proposition, 

\yN-XN\>lO~P'{\xN\ + \yN\) (22) 

follows. Condition can also be written as S > 2M10~P«e"^(^(^)-^). 

Consider now some precision m, the algorithm actually is working with. 
Assume for simplicity further that for the initial value X — Xq holds and 
assume without loss of generality S < ^uIp{xq). Then, first, yo = xq holds. 
Second, the above condition gives 2~^"^~^^^ > lO~Poe~^^^^^^~^^ since uIp{xq) < 
2-m+i|^^| < 2~'^~^^M holds. In other words, the above condition gives an 
upper bound m < pQ-\d{10)+N{\{x)—e)/ ln(2) — 1 on the needed precision m. 
Assume furthermore that m is big enough such that xn and ^at is computed 
with the demanded precision, that is \xn — xn\ < 10~P''|xAr| and \yN — yN\ < 
10^P°\yN\ holds. Using (122|) gives x^ 7^ ijN- But this is a contradiction 
since ^0 = ^^o- So the upper bound on m calculated above is still too small. 
Hence, m > po ■ Id(lO) + N{X{x) — e)/ ln(2) — 1 must hold. Since Condition 
( I2TI) also holds for any N' > N and the same po as well as for any p > po, 
(t{x,p) > (A(x) — £:)/ln(2) follows for all p > po. Computing a{x) finally 
gives the assertion. □ 

Furthermore, if X{x) = X{x) holds, then Corollary 13.31 gives (j{x) = 
max(0, A(x)) for the algorithm presented at the beginning of this sec- 
tion. Using the above proposition then leads to the following theorem. 

Theorem 3.3. Let {D, f) be as in Assumption \3. 1\ and x G QflZ^ given such 
that X{x) exists and X{x) = X{x) holds. Consider an algorithm computing an 
initial segment {xn)o<n<N of the orbit of x with relative error < 10~^. Then 
this algorithm has a loss of significance rate greater or equal to that of the 
algorithm specified by the recursion /[T^) and p^ . 
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4. Conclusions 

In this paper, two main issues are addressed. First it is shown that a 
mathematically rigorous treatment of the computability aspects of the itera- 
tion of a real function in terms of arbitrary-precision floating-point arithmetic 
including automated error analysis is straightforward. Also, this treatment 
is in a manner which is familiar to people working in the field of numerical 
analysis or scientific computing and also for theoretical computer scientists. 
Furthermore, the approach does not only allow answers concerning the exis- 
tence of an algorithm which meets the requirements of computability theory, 
but it also allows a treatment of its space complexity in form of the loss 
of significance rate (which is actually the lookahead in Computable Analy- 
sis) and optimality as discussed in the preceding section. As a consequence, 
the approach here enables a motivated reader the real implementation and 
supports a practical performance analysis. 

Second, the results show that the Lyapunov exponent, a central quantity 
in dynamical systems theory, also finds its way into complexity theory, a 
branch in theoretical computer science. In dynamical systems theory, the 
Lyapunov exponent describes the rate of divergence in the course of time of 
initially infinitesimal nearby states. For two states having a small but finite 
initial separation, the Lyapunov exponent has only relevance for short time 
scales [Si]. The reason is that due to the boundedness of the phase space, 
any two different orbits cannot separate arbitrarily far away. However, the 
loss of significance rate shows that the Lyapunov exponent has on long time 
scales not only an asymptotic significance but also a concrete practical one. 
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